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We obtain a new decomposition of the Riemann-Liouville operators of fractional in- 
tegration as a series involving derivatives (of integer order). The new formulas are valid 
for functions of class C n , n £ N, and allow us to develop suitable numerical approxima- 
tions with known estimations for the error. The usefulness of the obtained results, in 
solving fractional integral equations and fractional problems of the calculus of variations, 
is illustrated. 
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1 Introduction 

Let a ljx(t) := J a x(r)dT. If the equation 



Abstract 




(1) 



is true for the n-fold integral, n 6 N, 



then 



ir i x{t) = ji 



i 




Interchanging the order of integration gives 
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Since, by definition, ([I]) is true for n = 1, so it is also true for all n £ N by induction. The 
(left Riemann-Liouville) fractional integral of x{t) of order a > is then naturally defined, 
as an extension of ([T]), with the help of Euler's Gamma function V: 



The study of fractional integrals ((2J) is a two hundred years old subject that is part of a 
branch of mathematical analysis called Fractional Calculus [9j [T3"l 116] . Recently, due to 
its many applications in science and engineering, there has been an increase of interest in 
the study of fractional calculus Fractional integrals appear naturally in many different 
contexts, e.g., when dealing with fractional variational problems or fractional optimal control 
[U El El H21 [T3]. As is frequently observed, solving such equations analytically can be a 
difficult task, even impossible in some cases. One way to overcome the problem consists to 
apply numerical methods, e.g., using Riemann sums to approximate the fractional operators. 
We refer the reader to [H El [TOj [T7j and references therein. 

Here we obtain a simple and effective approximation for fractional integrals. The paper is 
organized as follows. First, in Section [21 we fix some notation by recalling the basic definitions 
of fractional calculus. In Section [3] we obtain a decomposition formula for the left and right 
fractional integrals of functions of class C n (Theorems 13.31 and 13. 4p . The error derived by 
these approximations is studied in Section[4l In Section[5]we consider several examples, where 
we determine the exact expression of the fractional integrals for some functions, and compare 
them with numerical approximations of different types. We end with Section[6]of applications, 
where we solve numerically, by means of the obtained approximations, an equation depending 
on a fractional integral; and a fractional problem of the calculus of variations. 

2 Preliminaries 

We fix notations by recalling the basic concepts (see, e.g., [9]). 

Definition 2.1. Let x(-) be an integrable function in [a,b] and a > 0. The left Riemann- 
Liouville fractional integral of order a is given by 




(2) 




t <E [a, b] 



while the right Riemann-Liouville fractional integral of order a is given by 




If a > 0, P > 0, and x E L p (a, b), 1 < p < oo, then 



a 



I?J?x{t) = J? +, x{t) and t L%4x{t) = t L^x(t) 



almost everywhere. The equalities hold for all t G [a, b] if in addition a + /3 > 1. 
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3 A decomposition for the fractional integral 



For analytical functions, we can rewrite a fractional integral as a series involving integer 
derivatives only. If x is analytic in [a, b], then 



1 oo 



(-l) k (t-a) k 
(n! )Z_, (fc + a)*! 



a;W(t) 



(3) 



for all t £ [a, b] (cf. Eq. (3.44) in [13]). From the numerical point of view, one considers finite 
sums and the following approximation: 



-l) fc (i-a) 



x {k \t). 



(4) 



One problem with formula ([3]) is the restricted class of functions where it is valid. In 
applications, this approach may not be suitable. The main aim of this paper is to present 
a new decomposition formula for functions of class C n . Before we give the result in its full 
extension, we explain the method for n = 3. To that purpose, let x E C 3 [a, b\. Using 
integration by parts three times, we deduce that 



Ta ^ (t-a) a . , (t - a) a+1 , . , (t - a) a + 2 



r(a + l)~ v " / ' r(a + 2) ~ v ~' ' T(a + 3) 
By the binomial formula, we can rewrite the fractional integral as 

al?x(t) = A- — ^-x(a) + ^— — ^x'(a) + ^— — x"(a) 



X » + FT^ [\t-rr +2 x^(r)dr. 
r (" + 3) A 



r(a + i) 



r(a + 2) 
(i-a) 



r(a + 3) 



a+2 



+ 



r(a + 3)'Sr(-a-2)p!(t-a)f A 



r(p - a - 2) 



(t - afx^i^dr. 



The rest of the procedure follows the same pattern: decompose the sum into a first term plus 
the others, and integrate by parts. Then we obtain 



rn , , (t-a) a . , (t-a) Q+1 „ , (t-a) Q+2 ... , 



r(a + l) 

(t-a) a+2 " 



+ 



r(a + 2) v y T(a + 3) 
r(p - a - 2) 



r(a + 2) S r (_ Q _i)(p_i)i( t _ ) P 



(t - a) a , , ft - a) Q+1 ,, . 



r(a + i) 

(i - a) a+2 



r(a + 2) 



" r(p-a-2) 
^ T(-a- 2)p! 

(r - a) p " V(tMt 
r(p - a - 2) 



p=2 



r(-a-l)(p-l)! 



+ 



+ 



*"(i) 



T(a + 3) 

( t _ a )°+2 °^ 
r(a + l) ^ 



r(p - a - 2) 
^ r(-a - 2)p\ 

Tjp- a- 2) 
+ 1 ) f- 2 T(-a)(p-2)\(t-a)P ./„ 



(r - a) p - 2 x'(r)dr 
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x(t) 



(t-a) a 
r(a + l) 

(t - a) a+2 
T(a + 3) 



p=3 



r(p- a- 2) 
r(-a)(p-2)! 



(t-a) Q+1 ,. . 
+ ^r, — -^x\t) 



T(a + 2) 



p=2 



r(p - a - 2) 
r(-a-l)(p-l)! 



+ 



, (t-ar +2 g. 



r(p-a-2) 

^ r( _ a _2) p ! 

r(p - a - 2) 



r(«) ^ T(-a + l)(p - 3)!(t - a)f J a 



(t - a) p " 3 x(r)(ir. 



Therefore, we can expand a Ifx{t) as 
al?ar(t) = A Q (a)(t - a) a x{t) + Ai(a)(t - a) a+ V(i) + A 2 (a)(i - a) a+ V'(t) 



where 



4) (a) 
Ai(a) 
^ 2 (a) 



1 



I> + 1) 
1 

" I> + 2) 
1 

~ r(a + 3) 
B(a,p) 



p=3 



p=2 

oo 



+ ^B(a,p)(t-ar+ 2 ~%(t), (5) 
p=3 



r(p - a - 2) 
r(-a)(p-2)! 

r(p - a - 2) 
r(-a-l)(p-l)! 

r(p - a - 2) 
r(-a - 2)p! 



p=l 
r(p - a - 2) 



r(a)r(l-a)(p-2)!' 



and 



f P (t)= [\p-2)(r-ay- 3 x(T)dr. 

J a 



(6) 



(7) 



Remark 3.1. Function V p given by ((?]) may 6e defined as the solution of the differential 
equation 

f ^(t) = (p-2)(t-a)P- 3 x(t) 
I W = 

/or p = 3, 4, . . . 

Remark 3.2. When a is not an integer, we may use Euler's reflection formula (cf. J^) 

r(a)r(i-a) 



7T 



sm to 



to simplify expression B(a,p) in ©. 

Following the same reasoning, we are able to deduce a general formula of decomposition 
for fractional integrals, depending on the order of smoothness of the test function. 
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Theorem 3.3. Let n £ N and x 6 C n [a,b\. Then 



n-l 



l I«x{t) = Y J M<x)(t 



{t) + Y j B(a,p)(t-aY +n - l -vV p (t), 



(8) 



i=0 



p=n 



where 



Ai(a) 
B(a,p) 



r(a + * + l) 



!+ E 



p = Tl — % 

r(p — a — n + 1) 
r(a)r(l -a)(p-n + 1)1 



r(p — a — n + 1) 
T(— a — — n + 1 + i)\ 



i = 0, 



,n— 1, 



(9) 



and 



^(t) = / (p-n + l)(r-a) p - n x(r)(iT, 



(10) 



p = n, n + 1, . . . 



A remark about the convergence of the series in Ai(a), for i G {0, . . . , n — 1}, is in order. 
Since 

T(p — a — n + 1) r(p — a — i) 



E 



E 



1 



r(— a — — n + 1 + i)l T(—a — i)p\ 

p=n—i v / \j / p = Q \ > L 

= 2 F 1 (-a-i,-,-,l) - 1, 



(11) 



where 2-^1 denotes the hypergeometric function, and because a + z > 0, we conclude that 
(|lip converges absolutely (cf. Theorem 2.1.2 in [S]). In fact, we may use Eq. (2.1.6) in [S] to 
conclude that 

ET(p — a — n + 1) 
f( 

p=n— i 



-a — i)(p — n + 1 + i)\ 



1. 



Therefore, the first n terms of our decomposition ([8]) vanish. However, because of numerical 
reasons, we do not follow this procedure here. Indeed, only finite sums of these coefficients are 
to be taken, and we obtain a better accuracy for the approximation taking them into account 
(see Figures 5(a) and 5(b)). More precisely, we consider finite sums up to order N, with 
N > n. Thus, our approximation will depend on two parameters: the order of the derivative 
n E N, and the number of terms taken in the sum, which is given by N. The left fractional 
integral is then approximated by 



n-l 



a I?x(t) « J^(a,JV)(i - a) a+l x^(t) + Y J B{a,p){t - a) 01 ^ 1 -^), (12) 



i=0 



p=n 



where 



r(a + i+l) 



N 



!+ E 



p=n—t 



T(p — a — n + 1) 
r(— a — i)(p — n + 1 + i) 



(13) 



and B(a,p) and V p (t) are given by © and (fT0|h respectively. 
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To measure the errors made by neglecting the remaining terms, observe that 

I ^> T(p — a — n + 1) 1 Y(p — a — i) 



TT E 



E 



T(a + i + l) V r(-a-i)(p-n + l + i)! r(a + i + l) ^ r(-a - i)pl 



r(a + i + l) 
-1 



' 2 F\{-a - i, -, -, 1) 

N-n+i+l 

El (p — a — zj 
. , . . , r(-a - i)p! 

p=0 



JV-n+l+i 



E 

p=0 



r(p — a — i) 
r(— a — i)p\ 



(14) 



Similarly, 



r(p — a — n + 1) 
r(a)r(l - a) (p-n + 1)! 



E 



r(a)r(i - a) 



7V-n+l 

E 

p=0 



p! 



(15) 



In Tables [T] and [2] we exemplify some values for (|14[) and (I15p . respectively, with a = 0.5 and 
for different values of N, n and i. Observe that the errors only depend on the values of N — n 
and i for (11411 , and on the value of N — n for (1151) . 



~"\^iV - n 
i 





1 


2 


3 


4 





-0.5642 


-0.4231 


-0.3526 


-0.3085 


-0.2777 


1 


0.09403 


0.04702 


0.02938 


0.02057 


0.01543 


2 


-0.01881 


-0.007052 


-0.003526 


-0.002057 


-0.001322 


3 


0.003358 


0.001007 


0.0004198 


0.0002099 


0.0001181 


4 


-0.0005224 


-0.0001306 


-0.00004664 


-0.00002041 


-0.00001020 


5 


7.124 x 10~ 5 


1.526 x 10" 5 


4.770 x 10~ d 


1.855 x ID-* 


8.347 x 10- y 



Table 1: Values of error (1141) for a = 0.5. 



N -n 





1 


2 


3 


4 




0.5642 


0.4231 


0.3526 


0.3085 


0.2777 



Table 2: Values of error (1151) for a = 0.5. 



Everything done so far is easily adapted to the right fractional integral. In particular, one 
has: 

Theorem 3.4. Let n G N and x £ C n [a,b]. Then 

71—1 OO 

t I%x{t) = £ Ma)(b - t) a+i x^(t) + J2 B(a,p)(b - t) a+n - l ^W p {t) , 



i=0 



p=n 



where 



Ai(a) 
B(a,p) 



(-If 



r(a + i + l) 



i+ E 



p=n— i 

(-l) n T(p-a-n + l) 
r(a)r(l -a)(p-n + l)V 



T(p — a — n + 1) 
. r(— a — i)(p — n + 1 + i) 



W p (t)= / (p-n + l)(b-T) p - n x( T )dT. 
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4 Error analysis 

In the previous section we deduced an approximation formula for the left fractional integral 
(Eq. (|12|) ). The order of magnitude of the coefficients that we ignore during this procedure 
is small for the examples that we have chosen (Tables Q] and [2]). The aim of this section is to 
obtain an estimation for the error, when considering sums up to order N. We proved that 



J?x{t) = — - 2 r rx(a) + • • • + — —x (n ij 



r(a + i)' 



T(a + n) 



(a) 



Via + n) 



a+n— 1 



+ (< - a) " + ":' /7l-^Y «W(r)dr. 



t — a 



Expanding up to order N the binomial, we get 



t — a 



t — a 



a+n— 1 



N 



ET(p - a- n+l) ( t - a\ p 
rn - ~ -■ ( ~ — - i + R ^ T ) 

p=0 



"(1 — a — n)p\ \t — a 



where 



Rn(t) 



E 1 \P ■ 
T(l 



r(p — a — n + l) / t — a 



(1 — a — n) p\ \t — a 



p=N+l 

Since r € [a,t], we easily deduce an upper bound for Rn(t): 



\Rn(t)\< Yl 

p=N+l 



T(p — a — n + 1) 



r(l — a — n)p 



a + n — 1 
p 



(■oo e (a+n— l) 2 +a+n— 1 
/ ; dp 



N 



P 



a+n 



- E 

p=N+l 
g (a+n-l) 2 +a+n-l 

(a + n - 1) jV"+-™-i 



£ E 

p=N+l 



Q (a+n-l) 2 +a+n-l 



„a+n 



Thus, we obtain an estimation for the error Et T {-)' 

U — a ^a+n g (a+n-l) 2 +a+n— 1 



\E tr (t)\ < L n 



r(a + n)(a + n - l)jV"+«-i ' 



where L n = max \x^ n \r) 

T€[a,t] 



5 Numerical examples 

In this section we exemplify the proposed approximation procedure with some examples. In 
each step, we evaluate the accuracy of our method, i.e., the error when substituting a It x by 
the approximation alfx. For that purpose, we take the distance given by 



E 



J?x{t) - J?x{t)\ dt 



Firstly, consider x\(t) = t 6 and x^it) = t with t G [0, 1]. Then 



0/°- 5 *!(t) 



r(4) 
r(4.5) 



r- 5 and oirx 2 {t) 



£(ii) 
r(n.5)' 



40.5 



7 



(cf. Property 2.1 in \9\). Let us consider Theorem 13.31 for n = 3, i.e., expansion ([5]) for 
different values of step N. For function x\, small values of N are enough (N = 3,4,5). For 



%2 we take N = 4,6,8. In Figures 1(a) and 1(b) we represent the graphs of the fractional 



integrals of x\ and X2 of order a = 0.5 together with different approximations. As expected, 
when N increases we obtain a better approximation for each fractional integral. 



0.7 1 

Analitic 

N=3, E=0.0071256 

0.6 N=4, E=0.0021377 

• N.5, E=0.0008907 




(a) oI? 5 (t 3 ) 



Analitic 

N=4, E=0. 076374 

N=6, E=0.018516 

• N=8, E=0. 00701 17 



(b) ol^it 1 



Figure 1: Analytic versus numerical approximation for a fixed n. 

Secondly, we apply our procedure to the transcendental functions x^{t) = & and x^{t) 
sin(t). Simple calculations give 



)/ t °- 5 x 3 (i) = Vtj2 



k=0 



T(k + 1.5) 



and lf- 5 x 4 (t) = Vt^2 



■1)H 



k+2k+l 



k=0 



r(2fe + 2.5) 



Figures 2(a) and 2(b) show the numerical results for each approximation, with n = 3. We see 



that for a small value of one already obtains a good approximation for each function. 



- Analitic 

N=3, E=0.0025506 

- N=4, E=0.000781 
N=5, E=0.00033089 




- Analitic 

N=3, E.0. 000861 14 

- N=4, E.0. 00025326 
N=5, E.0. 0001 0397 




(a) o/f-V) 



(b) o/ t °' 5 (sin(f)) 



Figure 2: Analytic versus numerical approximation for a fixed n. 



For analytical functions, we may apply the well-known formula @. In Figure [3] we show 
the results of approximating with Q, N = 1,2,3, for functions xs(t) and x&(t). We remark 
that, when we consider expansions up to the second derivative, i.e., the cases n = 3 as in (JS]) 
and expansion ((3J) with N = 2, we obtain a better accuracy using our approximation ((SJ) even 
for a small value of N. 



- Analitic 

N-1, E.0.090457 

- N.2, E=0.019758 
N.3, E=0.0035792 





(a) otf-V) 



(b) o/ t °' 5 (sin(f)) 



Figure 3: Numerical approximation using @ of previous literature. 

Another way to approximate fractional integrals is to fix N and consider several sizes for 
the decomposition, i.e., letting n to vary. Let us consider the two test functions xi(t) = t 3 and 
X2(t) = t 10 , with t £ [0, 1] as before. In both cases we consider the first three approximations 
of the fractional integral, i.e., for n = 1,2,3. For the first function we fix iV = 3, for the 
second one we choose N = 8. Figures 4(a) and |4(b)] show the numerical results. As expected, 
for a greater value of n the error decreases. 



- Analitic 

n=1, E=0. 024583 

- n.2, E.0. 012826 
n=3, E=0. 0071 256 




Analitic 

n-1, E=0.011389 

n.2, E=0.0078954 

• n=3, E=0.0070117 




(a) I? Ht 3 ) 



(b) ol^it 1 



Figure 4: Analytic versus numerical approximation for a fixed N. 
We mentioned before that although the terms Ai are all equal to zero, for i £ {0, . . . , n— 1}, 



we consider them in the decomposition formula. Indeed, after we truncate the sum, the error 
is lower. This is illustrated in Figures 5(a) and 5(b) , where we study the approximations for 
I$- 5 xi(t) and I?- 5 x 2 (t) with Ai^O and A» = 0. 




Figure 5: Comparison of approximation ([5]) and approximation with = 0. 



6 Applications 

In this section we show how the proposed approximations can be applied into different sub- 
jects. For that, we consider a fractional integral equation (Example 16. ip and a fractional 
variational problem in which the Lagrangian depends on the left Riemann-Liouville frac- 
tional integral (Example 16. 2p . The main idea is to rewrite the initial problem by replacing 
the fractional integrals by an expansion of type (EJ) or (|8j) , and thus getting a problem involving 
integer derivatives, which can be solved by standard techniques. 

Example 6.1 (Fractional integral equation). Consider the following fractional system: 

f oI^ 5 x(t) = 24 ^ 4 flfil 

\ x(0) = 0. { ' 

Since o^°' 5 i 3 ' 5 = r ^ 5 ) t a , the function t \-t t 3 5 is a solution to problem (fl"6j) . 

To provide a numerical method to solve such type of systems, we replace the fractional 
integral by approximations (J1J) and (|12p . for a suitable order. We remark that the order 
of approximation, N in @ and n in (|12p . are restricted by the number of given initial or 
boundary conditions. Since (|16p has one initial condition, in order to solve it numerically, we 
will consider the expansion for the fractional integral up to the first derivative, i.e., N = 1 in 
(jH) and n = 2 in (I12D . The order N in (]12p can be freely chosen. 

Applying approximation (J3|), with a = 0.5, we transform (|16p into the initial value problem 

f 1.1285t a5 x(t) - 0.3761t L5 x'(t) = ^i 4 , 
\x(0) = 0, 
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which is a first order ODE. The solution is shown in Figure 6(a) It reveals that the approxi 



mation remains close to the exact solution for a short time and diverges drastically afterwards. 
Since we have no extra information, we cannot increase the order of approximation to proceed. 

To use expansion ([8]), we rewrite the problem as a standard one, depending only on a 
derivative of first order. The approximated system that we must solve is 



' A)(0.5, N)t°- 5 x{t) + ^(0.5, N)t 1 - 5 x'(t) + J2p=2 B (°^P) tl - 5 ~ PV p( t ) 
V;(t) = (p-l)tP- 2 x(t), p = 2,3,...,N, 
x(0) = 0, 
I V p (0) = 0, p = 2,3,...,N, 



T(4.5) ,4 
24 L ' 



where Aq,A\ are given as in (|13p and B is given by Theorem 13.31 Here, by increasing N, we 
get better approximations to the fractional integral and we expect more accurate solutions 
to the original problem (|16p . For N = 2 and N = 3 we transform the resulting system of 
ordinary differential equations to a second and a third order differential equation respectively. 
Finally, we solve them using the Maple built in function dsolve. For example, for N = 2 the 
second order equation takes the form 

^"(*) = J V 2(t) + JM*) " 5.1542t 2 ' 5 
V(0) = 
V'(0) = x(0) = 

and the solution is x(t) = V^it) = 1.34t 3 5 . In Figure [6(b)| we compare the exact solution with 
numerical approximations for two values of N. 




(a) Approximation by ((4} 



Analitic 

N=2, E=0.12172 

• N=3, E.0. 




(b) Approximation by (JTJ 



Figure 6: Analytic versus numerical solution to problem (|16p . 



Example 6.2 (Fractional variational problem). Let a £ (0,1). Consider the problem 

J[ x (.)} = f ( I?x(t) - tf dt — > min, 
Jo 

T(a + 1.5) 



(17) 



x(0) = 0, x(l) 



r(i.5) 



n 



Problems of the calculus of variations of type (|17p . with a Lagrangian depending on fractional 
integrals, were first introduced in [3]. In this case the solution is rather obvious. Indeed, for 



x(t) 



r(1.5 + a) 



Vt 



(18) 



r(i.5) 

one has oIfx{t) = t and J[x(-)] = 0. Since functional J is non-negative, f|18f) is the global 
minimizer of (|17|) . 

Using (fT2|) . we approximate problem (fT7|) by 



M-)\ 



N 



A (a,N)t a x(t) + A 1 {a,N)t 1+a 3^{t) + ^B{a,p)t l+a ~ p V p {t) - t 

p=2 



dt 



mm, 



v;{t) = { P -\y- 2 x(t), P = 2,3 

V P (Q) = Q, p = 2,3,...,N, 
r(a + 1.5) 

r(i-5) ' 



)*** 5 3 



N, 



x (0) = 0, x(l) 



(19) 

This is a classical variational problem, constrained by a set of boundary conditions and 
ordinary differential equations. One way to solve such a problem is to reformulate it as an 
optimal control problem [15]. Let us introduce the control variable 

N 

u(t) = A {a,N)t a x{t) + A^NY^x'it) + ^B{a,p)t 1+a - p V p {t). 

p=2 

Then (|19|) becomes the classical optimal control problem 



J[x(-)] 



[u(t) - t] dt 



mm, 



N 



r'(t) = -AoA-h^xit) - ^A^ l B p r p V p {t) + A^ 1 t~ 1 ~ a u(t), 

p=2 

V p l (t) = (p-l)t p ~ 2 x(t), p = 2,3,...,N, 
V p (0) = 0, p = 2,3,...,N, 



x(0) = 0, x(l) 



T(q + 1.5) 

r(i.5) ' 



where A{ = Ai(a,N) and B p = B(a,p). For a = 0.5 and N = 2, application of the 
Hamiltonian system with multipliers Ai(i) and X 2 (t) |15| . gives the two point boundary value 
problem 



x'(t) = -AoA[ x t- l x{t) - AT/ 1 B 2 t- 2 V 2 {t)t 1 ~ a - i^- 2 t- 2 - 2Q Ai(t) + A^t" a , 



< x' 

V'(t) =x(t), 



X[(t) =A A^ 1 t- 1 X 1 (t)-X 2 (t), 
{ X' 2 (t) =A^ 1 B 2 t- 2 X 1 (t), 



(20) 



with boundary conditions 



s(0) = 0, 
V 2 (0) = 0, 



x(l) = l, 
A 2 (l) = 0. 
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In practice the value of n is chosen taking into account the number of boundary conditions. 
For problem (|17p . to avoid lack of boundary conditions, we take n = 2. On the other hand, the 
value of N is only restricted by the computational power and the efficiency of the numerical 
method chosen to solve the approximated problem. Figure [7] shows the solution to system 
(|20|) together with the exact solution to problem (|17|) and corresponding values of J. 



Analytic: J=0 

Approximation: N=2, J=4.4517e-016 




0.2 0.4 0.6 0.8 1 

t 



Figure 7: Analytic versus numerical solution to fractional variational problem (|17p . 
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